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During the last few years progress has been made on several fronts making it possible to revisit 
Cauchy-perturbative matching (CPM) in numerical relativity in a more robust and accurate way. 
This paper is the first in a series where we plan to analyze CPM in the light of these new results. 

One of the new developments is an imderstanding of how to impose constraint-preserving boimd- 
ary conditions (CPBC); though most of the related research has been driven by outer boimdaries, one 
can use them for matching interface boundaries as well. Another front is related to numerically stable 
evolutions using multiple patches, which in the context of CPM allows the matching to be performed 
on a spherical surface, thus avoiding interpolations between Cartesian and spherical grids. One way 
of achieving stability for such schemes of arbitrary high order is through the use of penalty techniques 
and discrete derivatives satisfying summation by parts (SBP). Recently, new, very efficient and high 
order accurate derivatives satisfying SBP and associated dissipation operators have been constructed. 

Here we start by testing all these techniques applied to CPM in a setting that is simple enough 
to study all the ingredients in great detail: Einstein's equations in spherical symmetry, describing a 
black hole coupled to a massless scalar field. We show that with the techniques described above, 
the errors introduced by Cauchy-perturbative matching are very small, and that very long term and 
accurate CPM evolutions can be achieved. Our tests include the accretion and ring-down phase of 
a Schwarzschild black hole with CPM, where we find that the discrete evolution introduces, with a 
low spatial resolution of Ar = M/ 10, an error of 0.3% after an evolution time of 1, 000, 000 M. For a 
black hole of solar mass, this corresponds to approximately 5 s, and is therefore at the lower end of 
timescales discussed e.g. in the coUapsar model of gamma-ray burst engines. 

PACS numbers: 04.25.Dm, 04.25.Nx, 04.70.Bw 



I. INTRODUCTION 

It is generally expected that the geometry of compact 
sources should resemble flat spacetime at large enough 
distances. This is true not only qualitatively, but through 
very precise f alloff conditions that are built into the for- 
mal definition of asymptotic flatness. Within this def- 
inition, the deviations from flat spacetime are well de- 
scribed (in the sense of the leading order behaviour of 
an expansion in powers of "1/r") by perturbations of 
the Schwarzschild spacetime 

Such perturbations can in turn be studied through the 
gauge invariant Regge-Wheeler and Zerilli (RWZ) for- 
malisms 1 2, 3]. These allow one to derive, after a spher- 
ical harmonic decomposition (that is, for each "{£, m)"), 
two master evolution equations for the truly gauge in- 
variant, linearized physical degrees of freedom. Due to 
the multipole decomposition, these equations involve 
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only one spatial coordinate (the radial one). The fact 
that they are one-dimensional implies that these mas- 
ter equations can be solved for very large computational 
domains with very modest computational resources. On 
the other hand, three-dimensional Cauchy codes are 
very demanding on their resource requirements. Even 
though mesh refinement can help in this respect, there 
is a limit to how much one can coarsen the grid in the 
asymptotic region; this limit is set by the resolution re- 
quired to reasonably represent wave propagation in the 
radiative zone. The use of a grid structure adapted to the 
physical geometry (possibly through multiple patches) 
can also help 's',^, but one still ends up imposing 
artificial (even if constraint-preserving) boundary con- 
ditions at the outer boundary. For example, one in gen- 
eral misses information about the geometry outside the 
domain [7J. 



Two approaches that at the same time provide wave 
extraction, physically motivated boundary conditions, 
and extend the computational domain to the radiative 
regime are Cauchy-characteristic f^, '9*1 and Cauchy- 
perturbative matching (CPM) |10, 11, 12]; this paper is 
concerned with the latter. The idea is to match at each 
timestep a fully non-linear Cauchy code to an outer one 
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solving, say, the RWZ equations . 

This paper is the first one in a series where we plan 
to revisit CPM in the light of some recent technical de- 
velopments — which we describe below — that should 
help in its implementation. Before discussing these de- 
velopments, we point out and summarize some features 
present in the original implementation of CPM which 
we hope to improve on: 

1. The non-linear Cauchy equations were solved on 
a Cartesian, cubic grid. On the other hand, the 
RWZ equations use a radial coordinate for the spa- 
tial dimension. Mixing Cartesian coordinates with 
spherical ones leads to the need for interpolation 
back and forth between both grids. Especially 
when using high order methods, this type of inter- 
polation might not only be complicated but also 
subtle: depending on how it is done it might in- 
troduce noise and sometimes it might even be a 
source of numerical instabilities. 

2. When injecting data from the perturbative mod- 
ule to the Cauchy code and vice versa bound- 
ary conditions were given to all modes, irrespec- 
tively of their propagation speed and without tak- 
ing into account the existence of constraint violat- 
ing boimdary modes. One would intuitively ex- 
pect a cleaner matching if boundary conditions are 
given according to the characteristic (propagation) 
speeds of the different modes, and even cleaner 
if constraint-preservation is automatically built in 
during the matching. 

3. Low order numerical schemes, which result in 
slow convergence, were used. 

In recent years there has been progress on several re- 
lated fronts that should in principle help in the imple- 
mentation of CPM. We describe these new results next^: 

1 . The first improvement is the ability to implement 
smooth (in particular, spherical) boundaries in 3D 
Cauchy evolutions |4, 5, 6]. One important advan- 
tage of this is the fact that the matching can be per- 
formed — to either a perturbative or a characteris- 
tic outer module — without the need for interpo- 
lation between spherical and Cartesian grids. In 



^ Even though including the angular momentum of the background 
is a high order correction in terms of powers of 1/r, one might, in 
principle, try to solve for perturbations of Kerr spacetime (as op- 
posed to Schwarzschild). 

^ There is actually another ingredient: the use of a generalized pertur- 
bative formalism that allows for any (spherically symmetric) slicing 
of the background Schwarzschild metric 1 13]. However, since such 
ingredient will not appear in the simplified model that we look at in 
this paper, we skip its discussion here. 



that way a possible source of noise can be elimi- 
nated. It is now understood how to match differ- 
ent domains using schemes of arbitrary high or- 
der while at the same time ensuring numerical sta- 
bility. One way of doing so is through the use of 
multiple patches (much in the same way multiple 
charts are used in differential geometry), penalty 
terms and difference operators satisfying summa- 
tion by parts |4] (more about this below). This is 
the approach we shall explore here in the context 
of CPM^. 

2. The second improvement is the construction 
of constraint-preserving boimdary conditions 
(CPBC). Several efforts have by now reported 
numerically stable (in the sense of convergent) 
implementations of such boundary conditions 
for the fully three-dimensional non-linear Ein- 
stein's equations Furthermore, there 
have been reports in the context of Cauchy- 
characteristic matching that significant improve- 
ments are obtained when this type of boundary 
conditions are used in the matching |16]. With 
this in mind, we will test their use in CPM. 

3. Lastly, new, accurate and efficient high order dif- 
ference operators satisfying SEP and associated 
dissipative operators have been constructed re- 
cently IT^ ^3, 19, 20]. As mentioned above, in 
conjunction with certain penalty interface treat- 
ment such operators guarantee numerical stability 
when "glueing" together different computational 
grids. We will test these operators in the context of 
CPM. 

We have incorporated these techniques, i.e., high- 
order summation-by-parts finite differencing and dis- 
sipation operators, multiple coordinate patches with 
penalty inter-patch constraint-preserving boundary 
conditions and Cauchy-perturbative matching, into 
a spherically symmetric numerical code evolving the 
Einstein-Christoffel form of the field equations |21|, 
minimally coupled to a Klein-Gordon field. Using this 
tool, we can test the performance of the numerical meth- 
ods in a non-trivial, but easily reproducible and com- 
putationally inexpensive setting, and gain experience 
for three-dimensional applications. The evolutions pre- 
sented here model black holes with excision in isola- 
tion, under dynamical slicings, and black holes accreting 
scalar field pulses, which are used as a scalar analogue 
of gravitational radiation. 

The plan of this paper is as follows. In section [H] 
we introduce the continuum system and the numeri- 
cal techniques we have used. Results are presented in 



^ Regardless of whether matching is present or not, the use of multi- 
ple coordinate patches has advantages when modelling black holes 
through excision of the singularity from the computational domain. 
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section |^ where a black hole is evolved successively 
from simple settings, i.e., single-patch, isolated. Killing- 
field adapted gauges, to more involved ones including 
Cauchy-perturbative matching and scalar pulse accre- 
tion. Finally, in section Hvl we draw conclusions and 
give an outlook to future work. 



II. EQUATIONS AND METHODS 

A. Evolution equations and constraint-preserving 
boundary conditions 

In this paper we use the Einstein-Christoffel (EC) sys- 
tem 1 21] in spherical symmetry. We follow the notation 
of Ref. 1 22]; in particular, the densitized lapse is denoted 
by a = Ng~^^^, and a = a.r^ sin{0) is introduced for con- 
venience. Here, g is the determinant of the 3-metric and 
N the lapse function, while the 4-metric is written as 

ds^ = -N^dt^ + grr{dr + ^dtf + r^gridO^ + sin^ Odcj)^) 

The vacuum part of the evolution equations in spher- 
ical symmetry for this formulation constitute a symmet- 
ric hyperbolic system of six first order differential equa- 
tions. The vacuum variables are the two metric and ex- 
trinsic curvature components 

grr, gTr Krr, Kj, 

where the extrinsic curvature is written as 



Mode 



Kij = K„dr^ 



r^Kjide^ 



sin^ dd<p'^) 



plus two auxiliary variables needed to make Einstein's 
equations a first order system. These extra variables are 
defined as 



f — off 
Jrrr — ^ 



Irrfrl 
gT 



rT 



S't 



gT 

r 



In addition, a massless Klein-Gordon field is mini- 
mally coupled to the geometry L^,22J. The scalar field 
equation 



g^'VaVfeY = 

is converted into a first order system by introduction of 
the variables 



n 



N 



(^Y'-Y), 



Throughout this paper the 'prime' and 'dot' represent 
partial derivatives with respect to r and t, respectively. 

Constraint preserving boundary conditions are im- 
posed by analyzing the characteristic modes of the main 
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U7 = n + <Pgn^/^ 
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r<2M 
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left 
left 
left 
left 
left 



r>2M 
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left 
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left 
right 
right 
right 
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Table I: Characteristic modes for Einstein-Christoffel system 
in spherical symmetry, and their direction of propagation for a 
Schwarzschild spacetime in Painleve-GuUstrand coordinates 
with respect to the vector field 3,-. In this gauge, all modes 
are outflow at the inner boundary, if it is located at r < 2M, 
while boundary conditions have to be applied to the incoming 
modes Mj, M2, M3, M4 and Mg at the outer boundary, assuming is 
is located at r > 2M. 



and constraint evolution systems, as discussed in 123 
These modes and their associated characteristic speeds 
are summarized in TableU] For illustration purposes, we 
also show the direction of propagation of each mode in 
the Schwarzschild spacetime in Painleve-Gullstrand co- 
ordinates 1 24, 25, 26]. 

From Table U we notice that for the Schwarzschild 
spacetime there are four ingoing and two outgoing grav- 
itational modes at the outer boundary, and therefore ex- 
pect the same count to hold for perturbations thereof. 
Boundary conditions for the incoming modes Mi, M2 ^rid 
M4 are fixed by the CPBC procedure. Thus, the only 
free incoming modes are 1/3, which represents a gauge 
mode and Ug, which represents a physical one (see 1 23| 
for more details). Boundary conditions do not need to 
be specified at the inner boundary if it is located inside 
the event horizon, because all modes are outflow then. 



B. Cauchy-perturbative matching 

Since there is no radiative degree of freedom in spher- 
ically symmetric spacetimes, we use the massless Klein- 
Gordon field as a scalar analogue of gravitational waves. 
To emulate the setup of three-dimensional Cauchy- 
perturbative matching as closely as possible, the scalar 
wave is evolved on a fixed Schwarzschild background 
in a "perturbative" patch defined for r > r,,,, while the 
fully non-linear Einstein's equations are evolved in the 
"Cauchy" patch, defined for r G [re, r,„], where and 
r„i denotes the excision radius and the matching radius, 
respectively. 

The fact that we are using CPBC allows us to perform 
a clean matching. From the analytical point of view our 
matching works in the following way: As mentioned 
above, after the CPBC procedure, only two free char- 
acteristic modes are entering the Cauchy computational 
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domain (at r = rm), denoted by M3 and wg. Since in a 
very precise sense M3 is a gauge mode, we are free to give 
boimdary conditions to it in a very simple way: we just 
set it to its initial value. Regarding u^, we use the "per- 
turbative" value of the same quantity coming from the 
perturtative domain as counterpart, and communicate 
these two modes (how this is done at the numerical level 
is explained below). Similarly, there is only one charac- 
teristic mode entering the perturbative domain, which 
is the linearized version of U7. We therefore communi- 
cate the non-linear and linear versions of that mode as 
well. 



C. Discrete techniques 

Given a well-posed initial-boundary value problem 
for Einstein's field equations, we construct a stable and 
accurate discrete system by using operators satisfying 
the SBP property. In short, a finite difference operator, 
D, satisfies SBP on a computational domain [a,b] dis- 
cretized using grid points i — 1, . . . ,n and a grid spacing 
hii 

{u,Dv) + {v,Du) = {uv)\l (1) 

holds for all grid functions u, v. Here the scalar product, 
E is defined in terms of its coefficients dij by 

n 

{u,v) J2 UiVjCTij. (2) 

In this paper we use the new, efficient, and accurate 
high order SBP difference operators and associated dis- 
sipation operators constructed in Ref. |20]. Thus, as 
mentioned, this paper also serves as an extra test of 
those new operators. 

SBP operators are standard centered finite difference 
operators in the interior of the domain, but the stencils 
are modified to yield lower order operators in a region 
close to the boimdaries (at the boundary itself the stencil 
is completely one sided). There are several types of SBP 
operators depending on the properties of the norm. The 
simplest are the diagonal norm operators. They have 
the advantage that SBP is guaranteed to hold in several 
dimensions by simply applying the ID operator along 
each direction and that numerical stability can be guar- 
anteed by discrete energy estimates in a wide range of 
cases. The main disadvantage is that the order of the 
operator at and close to the boundary is only half the 
interior order. We denote the SBP operators by the inte- 
rior and boundary order and consider here the diagonal 
operators D2_i, D4_2, Dg.s and D8_4. The second type 
is the restricted full norm operators, where the norm is 
diagonal at the boundary but has a non-diagonal block 
in the interior. The advantage of these operators is that 
the order at and close to the boundary is only one order 
lower than in the interior, while the disadvantage is that 



schemes based on these operators may be unstable with- 
out the use of dissipation. The restricted full operators 
we use here are 04.3 and 05.5. 

If the computational domain is split into several 
sub-domains ("patches"), the discrete representation re- 
quires a stable technique to communicate the solution 
at inter-patch boundaries. We make use of a penalty 
method |27], which adds a damping term to the right 
hand side of the evolution equation at the boundary 
point in a way which retains linear stability. The method 
has a free parameter, called S in Ref. |27], which deter- 
mines how much the difference between characteristic 
fields on either side of the inter-patch boundary is pe- 
nalized. Different values of S result in different amount 
of energy dissipation at the inter-patch boundary and 
can in principle be chosen so that no energy is dissipated 
(this is marginally stable). Usually the value of 5 is cho- 
sen such that some dissipation of energy occurs. With 
constant values of 5 the amount of dissipation decreases 
with resolution. 



D. Numerical code 

For the purposes of this paper, a one-dimensional 
code which supports constraint-preserving boundaries, 
multiple grid patches, and the use of the aforemen- 
tioned high order SBP derivative and dissipation opera- 
tors has been developed. In addition, the code is able to 
reproduce the (single grid and without CPM matching) 
second-order methods of Ref. |23] for comparison. We 
use the methods of lines, and the time integration is per- 
formed by a 4th order Runge-Kutta method. The grid 
patches that we consider here are not intersecting, but 
touching. This implies, that each grid function is double 
valued at the patch interface coordinate since the SBP 
derivative operators are one sided at the boundaries. To 
ensure consistency without compromising (linear) sta- 
bility, we make use of a penalty method as described 
above. Constraint-preserving boundary conditions re- 
quire the calculation of derivatives of certain grid func- 
tions at the outer boundary, which we also obtain by us- 
ing the SBP derivative operators. 

In a black hole setting, the computational domain next 
to the excision boundary tends to quickly amplify high 
frequency noise, which can not be represented accu- 
rately on the discrete grid. This is especially true for 
high order accurate derivative operators. Thus, high or- 
der simulations of black holes need a certain amount 
of numerical dissipation to be stable. This dissipation 
is here provided by the SBP dissipation operators con- 
structed in Ref. 120]. The free parameters of these oper- 
ators, namely the coefficient of the dissipation and the 
extent of the transition region (for non-diagonal opera- 
tors), are found by numerical experiment. 
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III. RESULTS 

The numerical experiments presented in this sec- 
tion are set up to systematically test the performance 
of the new techniques in several situations of increas- 
ing difficulty. We start with a series of tests evolv- 
ing a Schwarzschild black hole in Painleve-Gullstrand 
coordinates with either a single patch or two patches 
matched via the penalty method, and compare the 
performance of all SBP operators with the second or- 
der finite-differencing method presented in [23]. Next, 
to test more dynamical situations, a gauge or scalar 
field signal is injected in a constraint-preserving man- 
ner through the outer boundary and accreted onto the 
black hole. A robust stability test is then performed 
with noise on the incoming gauge mode M3, and, with 
Cauchy-perturbative matching, on the scalar field mode 
Ug. Finally, a series of high-precision tests involving all 
techniques are presented, in which a black hole accretes 
a scalar field injected through the outer boundary of the 
perturbative patch. These simulations also include a test 
of the long-term stability and accuracy after accretion 
and ring-down. 



A. Schwarzschild black hole in Painleve-Gullstrand 
coordinates 

In our first series of tests, a Schwarzschild black 
hole is evolved with high-order accurate SBP opera- 
tors, constraint-preserving boundary conditions and ex- 
cision. Cauchy-perturbative matching is not used in 
these tests. To fix the coordinate system, we make use 
of the horizon-penetrating Painleve-Gullstrand coordi- 
nates 1 25, 26], and we fix the coordinate fimctions a and 
/5 of the previous section to their exact values. 

For all tests, the inner boundary is located well in- 
side the event horizon (more precisely, it is located at 
Ye = IM), which implies that all modes are outflow. 
Therefore, no boundary conditions may be applied at 
the excision boimdary. The exact boundary location is 
not crucial as long as it is inside the apparent horizon, 
but this choice facilitates comparison with tZSll . Also, in 
dynamical situations the apparent horizon location may 
move significantly on the coordinate grid, and to ensure 
outflow conditions at the inner boundary some penetra- 
tion into the black hole is of advantage. To match the 
setup of 1 23), we set the outer boundary to r lOM. 
To ensure well-posedness of the continuum problem, 
boimdary conditions should be applied to the incoming 
modes Mi, U2, M3, M4, and ug. However, three of these 
modes, namely Ui, 112, and M4, can be fixed by the use 
of constraint-preserving boundary conditions, as dis- 
cussed in Section mi which leaves the freely specify able 
gauge mode M3 and the scalar field mode i/g- Since in 
these initial tests we are only interested in obtaining a 
stationary black hole solution, the initial scalar field is 
set to zero, and the (scalar field) characteristic mode Mg 



is penalized to zero as well. The incoming gauge mode 
M3 is penalized to the exact solution. 

An error function 5M can be defined by use of the 
Misner-Sharp mass function Qj 



M(r) ^ 




(3) 



where then, if the black hole mass is denoted by M, 
SM{r) = (M(r) — M) / M. Since the same error measure 
and continuum system is used in |23], we can compare 
the different discrete approaches directly. 



1. One grid patch 

The computational domain r G [1, 10] is represented 
by one coordinate patch, which is exactly the same setup 
as in ref. |23]. In Figure ^ we compare for coarse and 
high resolutions, Ar = M/8, M/64, the performance 
of the methods used in ref. [23], namely second order 
spatial derivatives with fourth order Kreiss-Oliger dis- 
sipation (which is set to zero near the boundaries) and 
a third order extrapolation at the boundaries, with the 
SBP derivative and dissipation operators D2-1, D4_2, 
D4_3, D6_3, Dg-s and D8_4. The figure shows the evo- 
lution of the L2 norm of the Misner-Sharp mass error 
over an evolution time of 10, DOOM. In all cases dis- 
played there is a linear growth in the error after some 
time. This is an artefact of the discrete representation 
of the constraint-preserving boundary conditions. We 
have also performed tests with maximally dissipative 
boimdary conditions: these yield a discrete equilibrium 
after some time, and thus allow for evolutions of un- 
limited time. However, since these boundary conditions 
are not correct for most systems of practical interest, we 
only make use of this result to point out the source of 
the linear growth of errors observed, which converges 
away with increasing resolution. 

As soon as the error gets close to 1, the code encoun- 
ters an instability, which, in this case, is associated with a 
migration of the excision boundary outside of the black 
hole, and consequently ill-posedness of the continuum 
problem. While this migration could be theoretically 
avoided by choosing horizon-fixing dynamical coordi- 
nate conditions, a solution with this magnitude of error 
is, in any case, not of practical use. 

In the present numerical code, the SBP operators are 
also used as one-sided derivatives for determining the 
constraint-preserving boundary conditions, which sug- 
gests that the operator D2-1, which is only first order at 
the boundaries, will yield less accurate outer boundary 
conditions than the third order method in [23]. Figure 
n] clearly demonstrates this fact. However, the opera- 
tors D5_3, D6_5 and Dg_4 are significantly more accu- 
rate than the results presented in ref. [23], and already so 
at the coarsest resolution. Furthermore, at Ar = M/64 
the SBP operator Dg-s induces a solution error of less 
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time [M] 

Figure 1: Time evolution of the relative error in the Misner- 
Sharp mass function when evolving a Schwarzschild black 
hole in Painleve-GuUstrand coordinates with one grid patch, 
for different discrete methods. Two resolutions are displayed, 
corresponding to Ar = M/8 (upper panel) and Ar = M/64 
(lower panel). The result from the method presented ref. 1 23] is 
denoted by "second order ", while new results are marked by 
the SBP derivative and dissipation operators used. The high- 
order operators Dg_5 and D8_4 display superior performance 
already at the lowest resolution. 



than 10 (that is, four orders of magnitude smaller than the 
corresponding errors when using the second order method of 
[23] with the same resolution) within 10, OOOM, which ap- 
pears sufficiently accurate for many practical purposes. 

The long-term evolution of a Schwarzschild black 
hole with the operators Dg-s and is displayed in 
Figure |3 The linear growth of errors dominates the 
solution at late times, but since this error significantly 
decreases with resolution, long evolution times can be 
obtained even for moderate radial grid spacings. This 
is naturally an interesting feature for simulations with 
three-dimensional spatial grids, where computational 
resources are still a viable concern. 



D^_3, M/8 




time [M] 

Figure 2: Evolution of a Schwarzschild black hole for 
100, OOOM. The axes show the quantities described in Figure^ 
It is clear that even with low resolutions of Ar = M/8 and 
M/16, the operators D5_5 and Dg_4 are able to evolve the 
black hole in a stable manner for a significant time. 



2. Tivo grid patches 

As dicussed in the introduction, the use of multi- 
ple coordinate patches has advantages when modelling 
black holes. To implement a stable interface bound- 
ary condition, the penalty method is used to ensure lin- 
ear stability. Here we first investigate the performance 
of the SBP operators coupled to an inter-patch penalty 
boundary method by evolving a black hole spacetime 
covered by two non-intersecting spherical shells, the 
first one from r = IM to r = 5.5 M, and the second 
one from r = 5.5M to r = lOM. In order to provide 
an intermediate test towards the CPM tests below, we 
do a non-linear matching, communicating all charac- 
teristic modes (that is, without imposing for the mo- 
ment constraint-preserving boundary conditions at the 
matching interfaces). 

The free parameter of the penalty boundary condition 
S introduced in section lll Cl is set to the dissipative value 
0. Only the operators Dg.s and D8_4 are used for com- 
parison to the results from the previous section. 

In Figure |3] the performance of the multi-patch sys- 
tem is compared to the uni-patch results from the pre- 
vious section. As expected, the use of one-sided deriva- 
tives at the inter-patch boundary reduces the total level 
of accuracy, but in a very small amount; furthermore, 
the system is still stable and convergent. The time of 
the onset of the linear growth observed in all evolutions 
varies between the grid setups and choices of discrete 
operator. Figure HI shows the 3-metric component grri^") 
at the times t ^ and t = 10, OOOM. The region around 
the inter-patch interface at r = 5.5M is shown in the in- 
set, which demonstrates that the penalty method intro- 
duces no strong visible artifacts in this part of the solu- 
tion. This observation also holds for the other solution 
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Figure 3: Comparison of uni-patch and multi-patch evolutions 
of a Schwarzschild black hole in Painleve-GuUstrand coordi- 
nates. The graphs denoted by "one patch" and "second or- 
der" are those from Figure^ while the corresponding graphs 
for "two patches" cover the computational domain with two 
non-intersecting spherical shells, the first one from r = IM to 
r = 5.5M, and the second one from r — 5.5M to r = lOM. 
The one-sided derivatives at the interface boundary introduce 
a very small loss of accuracy. In the upper and lower panels 
the resolution is Ar = M/8, M/64, respectively. For the late 
time behaviour of D(,_^ and D8_4 please also cf. Figure|3 



fiinctions. 



B. Gauge wave on a Schwarzschild background 

The next series of tests focuses on a dynamical sit- 
uation, namely the evolution of a Schwarzschild black 
hole in non-stationary coordinates. For this purpose, 
the initial data is set to a Schwarzschild black hole in 
Painleve-Gullstrand coordinates as in section lllll A> , as 
is the lapse and shift function at all times, but the in- 
coming gauge mode M3 at the outer boundary is set to a 
Gaussian pulse of the form 



t = o 

t= 10,000 m 




r[M] 



"3(f) 



,PG 



{1 + Ae- 



-to)V^^ 



(4) 



r [M] 



Figure 4: Evolution of metric function grr for a black hole in 
Painleve-Gullstrand coordinates, with a resolution of Ar = 
M/64, two grid patches with an interface at r = 5.5M and 
using the SBP operator Dg_5. The two graphs show the metric 
function at t = (where grr(r) = 1) and at t = 10,000M. The 
inset shows the region around the interface between the grid 
patches. 



Here, 1(3^ is the exact gauge mode from the stationary 
solution. As in ref. |23], we impose a strong pulse with 
A = 1, to = 5M and a = 2M. Since the solution is 
now not adapted to the asymptotically timelike Killing 
field, the SBP operators and multi-patch techniques can 
be tested on a solution with wave propagation without 
compromising the use of the error measure | |i5M| |2- To 
facilitate comparison with ref. |23], the outer boundary 
is located at r = 30M in these tests. 

Figure |5] shows results from the gauge pulse problem 
on a single grid patch and two grid patches, here with 
an inter-patch boundary at r = 15.5M. While in the 
stationary case the inter-patch boundary method only 
had to deal with small numerically introduced differ- 
ences between the values of the geometrical quantities at 
the interface, the non-stationary case introduces a large 
pulse travelling over the boundary, and is thus a much 
more severe test for accuracy and stability of the penalty 
method. The solution error is dominated by the ability 
of the discrete method to represent the propagation and 
accretion of the gauge pulse, and by possible artefacts 
introduced by the inter-patch boundary. 

Judging from Figure |5l the high-order operators are 
stable and significantly more accurate than a second 
order method also in a dynamical situation, and even 
when using multiple matched domains. 



C. Accretion of a scalar field pulse 

Since the outer boundary has two free incoming 
modes, it is possible to inject a scalar fi eld pu lse in a 
way similar to the gauge pulse of section Jill B> . In con- 
trast to the gauge pulse, however, this system will result 
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Figure 5: Comparison of uni-patch and multi-patch evolutions 
of a gauge wave travelling on a Schwarzschild background. 
The graphs denoted by "second order " are obtained with the 
methods in ref. |23], while the corresponding graphs for "one 
patch" and "two patches" cover the computational domain 
with either one or two non-intersecting spherical shells, the 
first one from r = IM to r = 15. 5M, and the second one from 
r = 15. 5M to r = 30M. The one-sided derivatives at the in- 
terface boundary introduce a small loss of accuracy, but the 
system is still stable. The upper and lower panels correspond 
to Ar = M/8,M/64, respectively. 



in an increase of mass of the black hole, which also im- 
plies that the Misner-Sharp mass cannot be used as a 
measure of the errors anymore. A possible choice for a 
gauge field source with compact support is 







t < t, 



"8(0 







t > tp 



To facilitate comparisons with ref. [23] we use an am- 
plitude A = 7.2, and tj = OM, tp = lOM and set the 
computational domain to be r G [1,50]M. 

For resolutions Ar = M/20 and Ar = M/40, the time 
evolution of the apparent horizon is shown in Figure |6] 
The scalar pulse leads to a significant increase in the 



o ■ 9 ■ 0. O I 



Q-o M/20 
■ ■ M/40 



100 
time [M] 



Figure 6: Evolution of the apparent horizon mass for the ac- 
cretion of a strong scalar pulse to a Schwarzschild black hole. 
Shown are plots for two resolutions, Ar — M/20 and Ar = 
M/40, using the SBP operator D^_^. The large scalar field am- 
plitude leads to a significant increase in the black hole mass. 
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Figure 7: L2 norm of the Hamiltonian constraint over time for 
the accretion of a strong scalar field pulse to a Schwarzschild 
black hole, with resolutions Ar = M/20, M/40 (upper and 
lower panels, respectively). The graph denoted by "second 
order" is obtained with the method presented in |23], and the 
D5_5 and D8_4 are obtained using the corresponding SBP op- 
erators. 
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time [M] 

Figure 8: As Figure |7| but evolved for 10,000 M with Ar = 
M /20 to demonstrate the long-term behaviour after accretion 
of the pulse. 



black hole mass by a factor of w 2.7 after the pulse is 
inside the black hole. Larger amplitudes are not obtain- 
able with the simple gauge prescription used here, but a 
horizon-freezing gauge condition could improve on this 
result. As a replacement for the Misner-Sharp error mea- 
sure, we plot the L2 norm of the Hamiltonian constraint 
over time in Figure |7| It is apparent that the high-order 
operators are again stable and more accurate than the 
second order operator. The graphs indicate a growth of 
the constraint near f = 200M, but a long-term evolution 
with Ar = M/20 shown in Figure |H] demonstrates that 
the system settles down to stability after the accretion. 



D. Robust stability test with gauge noise 

The term robust stability test 12^ typically refers to the 
discrete stability of a numerical system in response to 
random perturbations. In this case, we will use the same 
system as in section JIII A2> , but impose random noise 
on the incoming gauge mode M3 with a certain ampli- 
tude. To test the discrete stability of the evolution sys- 
tem, we chose a large range of amplitudes from 10~* to 
0.3. Random perturbations of the latter amplitude is sig- 
nificant for a non-linear system*. 

For this multi-patch test, results in the mass error for a 
resolution Ar = M/8 are shown in Figure|9] It is appar- 
ent that strong random noise induces a stronger growth 
in the solution error. However, this growth is still lin- 
ear. As in all black hole evolutions in section lllll A> , the 



* Beyond this amplitude the inner boundary tends to become par- 
tially inflow by moving the apparent horizon beyond the compu- 
tational domain. More sophisticated gauge or inner boundary con- 
dition could alleviate this, but since we are interested here in a proof 
of principle, a simple system is preferred. 







1 




1 


' / 






— A = 














■ ■ A = 10"* 










; 




— ■ A= 10""^ 












- 


« A = 10"^ 














A= 10"' 














_1 

>^ A = 3 X 10 












r 


























: 














































■ ■ e ■ ■ — e- ■ ■ 0. 


■■■■<>■ 














1 











■4 1 , , I , , , , I 

10 100 1000 iOOOO 



time [M] 

Figure 9: Results of a robust stability test for different random 
noise amplitudes. The system is a Schwarzschild black hole 
in Painleve-GuUstrand coordinates, and the computational do- 
main r e [1, 10] M is covered by two patches with a boundary 
at r = 5.5M and a resolution of M/8. Random noise is super- 
imposed on the ingoing gauge mode U3, with an amplitude 
denoted by A. The graphs show the mass error with time for 
different random noise amplitudes, obtained with the SBP op- 
erator D5_5. 




10 100 1000 10000 

time [M] 

Figure 10: Like Figure|9| but for the highest random noise am- 
plitude 0.3 and different resolutions. 



system encoimters a numerical instability as the solu- 
tion error approaches 1, but this is not a consequence 
of the random noise, but of the inner boundary becom- 
ing partially inflow due to a coordinate motion of the 
apparent horizon. Also, with increasing resolution, the 
growth rate of the error does not increase, as shown in 
Figure^! We conclude that this high-order evolution 
system is discretely stable against strong random per- 
turbations. 
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time [M] time [M] 



Figure 11: Robust stability test with Cauchy-perturbative 
matching. The system is a dynamically evolved Schwarzschild 
black hole in Painleve-GuUstrand coordinates matched to a 
perturbative module at r = 5.5M as described in the intro- 
duction. Random noise is imposed via the incoming scalar 
field mode at the outer boimdary. Plotted is the L2 norm 
of the Hamiltonian constraint over time for different noise 
amplitudes. All evolutions were done with a resolution of 
Ar = M/8 and the SBP operator D(,_5. 



Figure 12: Like Figure ^3 but for the highest random noise 
amplitude and different resolutions. 



by-parts operators. 



F. Cauchy-perturbative matching: Accretion of a 
"gravitational wave" and long-term evolution 



E. Cauchy-perturbative matching: robust stability test 
with scalar field noise 

We now test the stability of the system with Cauchy- 
perturbative matching against random perturbations in 
the scalar field. To this end, the comp utational domain is 
again subdivided as in section Jill D> , but the right patch 
evolves the scalar field on a fixed Painleve-Gullstrand 
background as explained in the introduction. The inter- 
patch boundary is thus matching the Cauchy patch to 
a perturbative one, and we test the stability of the sys- 
tem against random perturbations by imposing random 
noise on the incoming scalar field mode on the outer 
boundary of the perturbative patch. 

Since the mass error is not available for a system ac- 
creting a scalar field, the L2 norm of the Hamiltonian 
constraint is used again in Figure El No exponential 
growth can be observed in the Hamiltonian constraint 
violation. The same is true when increasing the reso- 
lutions, as in Figure El which also deserves some ad- 
ditional comments: The robust stability test does not 
lead to a converging sequence of solutions if the ran- 
dom noise amplitude is not diminished with resolution. 
However, the purpose of these tests is to excite any un- 
stable high freqcency modes present in the numerical 
system. The absence of any mode growing with increas- 
ing resolution shows that the system with a Cauchy- 
perturbative matching interface is stable even against 
strong random noise injected into the system. This is a 
promising result for any effort to do three-dimensional 
matching between Cauchy modules and perturbative 
ones using multiple patches and high-order summation- 



Finally, using the massless Klein-Gordon field as a 
scalar analogue of gravitational waves in spherical sym- 
metry, we model the accretion of a gravitational wave 
packet across a Cauchy-perturbative matching bound- 
ary. This test is an extension of the single-patch scalar 
field accretion of section <IIICL and makes use of all in- 
gredients presented in this paper for a stable and accu- 
rate evolution of black holes with Cauchy-perturbative 
matching. 

Since Cauchy-perturbative matching assumes the 
gravitational wave to be a small perturbation of a 
fixed background in the wave zone, the amplitude of 
the wave packet that we inject through the outermost 
boundary is chosen to be A = 0.01. Similarly to sec- 
tion llll CI we describe the packet by the function 

{0 t<ti 
f > 

where for the number of half waves in the pulse we set 
n = 100. We inject the packet from tj = to if = 
lOOM. The plots in Figure E| display the the evolution 
of the grid function O, and specifically the behaviour 
of the fimction around the Cauchy-perturbative match- 
ing interface, which is at r = 25.5M. The correspond- 
ing increase in apparent horizon mass is shown in Fig- 
ure El The evolution of the Hamiltonian constraint vi- 
olation using the SBP operator Dg-s and different res- 
olutions is shown m Figure El ^ is apparent that 
with the techniques used not only is the discrete sys- 
tem stable and accurate, but also the amount of non-linear 
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Figure 14: Accretion of a scalar wave packet across a Cauchy- 
peturbative matching interface. This plot shows the apparent 
horizon mass over time for evolutions with different resolu- 
tions and the SBP operator Dg-s. 
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Figure 15: Accretion of a scalar wave packet across a Cauchy- 
peturbative matching interface. This plot shows the L2 norm 
of the Hamiltonian constraint for different resolutions, using 
the SBP operator D(,_^. The non-linear constraint violations in- 
troduced at the continuum by the matching are small enough that 
they cannot he detected in these very accurate simulations. Please 
note, for comparison with Figure|5| that the amplitude of the 
Klein-Gordon signal is smaller compared to section IIII Ct . 



Figure 13: Accretion of a scalar wave packet across a Cauchy- 
peturbative matching interface, as a scalar analog for gravita- 
tional wave accretion in three-dimensional simulations. The 
packet consists of 50 waves injected from t = to t = lOOM, 
as described in the text. Here, the grid function O is plotted 
over the radial coordinate at f = 30M, 65M, llOM (from top 
to bottom), for the resolution Ar = M/20 and the SBP opera- 
tor D5_5. The inset shows the behaviour of the grid function 
around the matching interface, which is at r = 25.5M. Note 
that even though the grid function is in principle two-valued 
on the interface, the penalties in conjunction with high order 
operators only lead to a very small mismatch. 



constraint violations introduced at the continuum by the 
Cauchy-perturbative matching are very small, in Figure [THI 
they must actually be smaller than 10~^. 

The advantages of using high-order methods is made 
evident in Figures 1161 1171 1181 and Q^l Ir* these plots, 
the performance of the SBP operator Dg-s, which is 
sixth order in the interior and fifth order at the bound- 
aries, is compared to that of the operator 04.3, which 
is fourth order in the interior and third order at the 
boimdaries, for different choices of resolution. Al- 
though both operators show convergence, for a mass 
increase of about 10~^, the operator 04.3 is unable to 
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Figure 16: Accretion of a scalar wave packet across a Cauchy- 
peturbative matching interface. To demonstrate the advantage 
of using high-order methods, {Mah — 1) is shown for evolu- 
tions obtained with the SBP operators D4_3 Dg-s, with res- 
olution Ar = M/10. The loss of mass after accretion of the 
wave packet with compact support in t e [0, 100] M is a nu- 
merical artefact, which converges away with resolution. The 
inset shows that the evolution obtained with the operator D^_^ 
is not unstable, but only significantly less accurate. 
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Figure 18: Like Figure^] but for a resolution of Ar = M/40. 
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Figure 19: Like Figure^] but for a resolution of Ar = M/80. 



Figure 17: Like figure^] but for a resolution of Ar = M/20. 



reproduce the correct behaviour with reasonable grid 
resolutions. We consider this specifically important for 
three-dimensional simulations, where the necessary re- 
sources scale with if n denotes the number of grid 
points in each direction. Thus, for all simulations requir- 
ing a certain amount of precision, high-order operators 
are an essential requirement. 

The long-term evolution of a Schwarzschild 
black hole accreting a wave packet over a Cauchy- 
perturbative matching interface and settling down to 
equilibrium is shown in Figure |20| The black hole is 
evolved for 1,000, OOOM with the lowest resolution 
Ar = M/10 and the SBP operator 05.5. While an 
evolution of this length might appear to be of only 



2x10 4x10 6x10 8x10 

time [M] 



Figure 20: Long-term stable evolution of a Schwarzschild 
black hole after accretion of a scalar wave packet with Cauchy- 
perturbative matching. The SBP operator D5_5 is used with a 
resolution of Ar — M/10. Plotted are the apparent horizon 
mass and the Hamiltonian constraint over time. The apparent 
horizon mass indicates that the discrete evolution introduces 
a relative error of about 0.3% after 1, 000, OOOM. 
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technical interest, we note that modelling phenomena 
like hypernovae and collapsars in general relativity will 
require the stable evolution of a black hole for at least 
several seconds, which is the lower end of timescales 
associated with the coUapsar model of gamma-ray 
burst engines Jz^ . For a stellar mass black hole, 
M = Mq « 5}LS, that is Is w 200, OOOM0. 

IV. CONCLUSIONS AND OUTLOOK 

To obtain long-term evolutions of compact astro- 
physical systems in three spatial dimensions, advanced 
numerical techniques are preferable in that they may 
improve stability and accuracy of the associated dis- 
crete model system. While high accuracy enables effi- 
cient use of the available computational resources, well- 
posedness of the continuum model and numerical sta- 
bility are requirements which can not be met by in- 
creasing computational power. A number of techniques 
has been suggested to address these issues |27]: Mul- 
tiple coordinate patches, typically adapted to approxi- 
mate symmetries of certain solution domains, combined 
with high-order operators are expected to increase the 
accuracy of any model of a stellar system. Cauchy- 
perturbative matching provides an efficient way to ac- 
curately model the propagation of gravitational waves 
to a distant observer, and to yield physical boundary 
conditions on incoming modes of the Cauchy evolution. 
Constraint-preserving boundary conditions isolate the 
incoming modes on the constraint hypersurface, and, fi- 
nally, for evolving black holes, an excision boundary is 
desirable to concentrate on the behaviour of the exter- 
nal spacetime. Only recently the consideration of the 
well-posedness of the differential system and the appli- 
cation of theorems on discrete stability of the numerical 
system have provided hints as how to address the out- 
standing issues. In this paper, we have applied all these 
techniques to a model system: a spherically symmetric 
black hole coupled to a massless Klein-Gordon field. 

We find that the use of a first-order hyperbolic for- 
mulation of Einstein's field equations, combined with 
high-order derivative and dissipation operators with 
the summation-by-parts property, penalized inter-patch 
boundary conditions and constraint-preserving outer 
boundary conditions leads to a stable and accurate dis- 



crete model. Specifically, isolated Schwarzschild black 
holes in coordinates adapted to the Killing fields, and 
in coordinates on which a gauge wave is imposed, and 
Schwarzschild black holes accreting scalar wave pulses 
were taken as typical model systems involving excision. 
The results show that the introduction of several coor- 
dinate patches and of a Cauchy-perturbative matching 
interface does not introduce significant artefacts or in- 
stabilities. Rather, the high-order methods allow the ac- 
curate long-term evolution of accreting black holes with 
excision and Cauchy-perturbative matching in reason- 
able resolutions. As an example, we have presented the 
evolution of such a system with the high-order SBP op- 
erator D6_5, which, at a resolution of Ax = M/10, intro- 
duced an error of only 0.3% after an evolution time of 
1,000, DOOM. 

Most system of interest in general relativistic as- 
trophysics will necessarily require the use of three- 
dimensional codes. Results from a one-dimensional 
study are useful in that they (i) allow to gain expe- 
rience in a clean but non-trivial physical system, (ii) 
can be easily reproduced without the need to imple- 
ment three-dimensional codes with multiple coordinate 
patches and (iii) allow to isolate sources of difficulty 
in the three-dimensional setting more easily. With the 
promising results from this study, we will, as a next step, 
apply these techniques to a three-dimensional, general 
relativistic setting. 
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